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ABSTRACT 



We present a combined analysis of previously published high-precision radial velocities and astrometry for the GJ 876 planetary 
system using a self-consistent model that accounts for the planet-planet interactions. Assuming the three planets so far identified in 
the system are coplanar, we find that including the astrometry in the analysis does not result in a best-fit inclination significantly 
different than that found by Rivera and collaborators from analyzing the radial velocities alone. In this unique case, the planet-planet 
interactions are of such significance that the radial velocity data set is more sensitive to the inclination of the system through the 
dependence of the interactions on the true masses of the two gas giant planets in the system (planets b and c). The astrometry does 
allow determination of the absolute orbital inclination (i.e. distinguishing between i and 180 - i) and longitude of the ascending 
node for planet b, which allows us to quantify the mutual inclination angle between its orbit and planet c's orbit when combined 
with the dynamical considerations. We find that the planets have a mutual inclination <l> bc = 5.0° + '"' . This result constitutes the first 
determination of the degree of coplanarity in an exoplanetary system around a normal star. That we find the two planets' orbits are 
nearly coplanar, like the orbits of the Solar System planets, indicates that the planets likely formed in a circumstellar disk, and that 
their subsequent dynamical evolution into a 2:1 mean motion resonance only led to excitation of a small mutual inclination. This 
investigation demonstrates how the degree of coplanarity for other exoplanetary systems could also be established using data obtained 
from existing facilities. 

Key words, stars: individual: GJ 876 - planetary systems - astrometry - methods: data analysis 



1. Introduction 

A planetary system's "architecture" consists of the census, 
masses, and orbits of the objects in the system. The formation 
and evolutionary history of a planetary system are encoded in 
these characteristics. The architecture of the Solar System is 
broadly consistent with the theory of planet formation in a cir- 
cumstellar disk, and its deviations from the expected pattern are 
interpreted as evidence of evolutionary processes. It is impor- 
tant to determine the architecture of exoplanetary systems so as 
to have broader constraints on the formation and evolution of 
planetary systems than can be obtained from study of the Solar 
System alone. 

The GJ 876 planetary system is one exoplanetary system that 
has been the focus of much attention to ascertain its architecture. 
GJ 87 6 itself is an M4 dwarf star with an es sentially solar metal- 
licity dBonfils et alj|2005t iBean et al.ll2006l) that is at least older 
than 1 Gyr as evid enced by its low activity an d slow rotation 
( Marcv et alJll998b . It is at a distance of 4.7 pc (Perrvma n et alj 
119971) and appears to be a singleton as no stellar companions 
have been reported. The star does not appear to also harbor a sub - 
stantial debris disk dTrilling et alJl2000HShankland et alj|2 008). 



* Based on observations made with the NASA/ESA Hubble Space 
Telescope, obtained at the Space Telescope Science Institute, which 
is operated by the Association of Universities for Research in 
Astronomy, Inc., under NASA contract NAS 5-26555 (programs GO- 
8102, 8775, and 9233); and on observations obtained at the W. M. Keck 
Observatory, which is operated jointly by the University of California 
and the California Institute of Technology. 



A gas giant planet was found to orbit GJ 876 b y two groups 
independently in 1998 usin g Doppler spectroscopy ( Marc v et al.l 
[T99alDelfosse et al.lfl998h . This was the first convincing detec- 
tion of a planet around an M dwarf, and GJ 876 still remains one 
of the few know n planet host i ng sta rs of this type. Subsequent 
observations by iMarcv et~ai1 d2001l) revealed that the system 
contains a second, lower-mass gas giant that is in a 2:1 mean 
motion resonance with the first detected planet. 

Soo n after the discovery of t he sec ond p lanet in the GJ 876 
system , lLaughlin & Chambers! d200lb and iRivera & Lissauerl 
(2001) pointed out that the two identified planets are experienc- 
ing mutual interactions on an unprecedentedly short timescale. 
They found that the resulting perturbations are of such signif- 
icance that a model based on Keplerian orbits was not accu- 
rate enough to match the radial velocity data, and instead di- 
rect integration of the equations of motion (i.e. Newtonian or- 
bits) for the three body configuration was needed. Those au- 
thors also found that, although they introduced significant ad- 
ditional complexity in the modeling of the observational data, 
the occurrence of short-timescale interactions could allow in- 
ference of the true masses of the planets from radial velocity 
data alone. This is in contrast to the usual situation where only 
planets' minimum masses are determinable from radial veloc- 
ity data owing the orbital inclination degeneracy. The key to the 
additional insight is that the non-Keplerian perturbations are de- 
pendent on the true masses of the interacting bodies. Therefore, 
if these perturbations could be characterized well enough with 
radial v elocity data then the planet mass es could be con- 
strained. Laug hlin & Chambers! (|2001) and IRivera & Lissauerl 
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(2001) both concluded that tight limits on the planet masses 
could not be set given the data available at the time, but that fu- 
ture analysis of the continuing Doppler measurements for GJ 876 
wo uld possibly give bette r results. 

iBenedict et al.l (12002) carried out astrometric observations of 
the GJ 876 system using the Fine Guidance Sensor (FGS) instru- 
ment on the Hubble Space Telescope (HST) beginning around 
the time the first planet was announced and continuing for 2.5 
years. Analysis of these data revealed a residual perturbation 
with semimajor axis of 0.3 +0.1 mas in phase with the orbit 
of planet b expected from modeling radial velocity data. This 
was the first definitive astrometric detection of an exoplanet, 
and it is still one out of only a few such successful detections. 
Modeling the astrometry together with the radial velocities avail- 
able yielded an estimate of the orbital inclination of planet b 
(ib = 84° + 6°), and thus the planet's true mass after assum- 
ing a mass for GJ 87 6 (nib = 1.9 ± 0.3 Mj up ). In their analysis, 
Bene dict et al.l d2002l) used the standard Keplerian rather than 
Newtonian orbital calculations for modeling the radial velocity 
and astrometry data. 

Continuing the trend of exopl anet firsts and rarities in the 
GJ 876 system, iRivera et al.l (120051) proposed the existence of a 
third, very low-mass planet (mj sin id — 5.9 + 0.5M ffi ) based on 
a relatively extensive set of new high-precision radial velocities. 
At the time, this was possibly (allowing for the inclination ambi- 
guity) the lowest mass planet yet found around a main sequence 
star, and it is still one of only a few known planets with a mass 
potentially in the "Super-Earth" regime (i.e. 1 % m % 10M ffi ). 
A comprehensive photometric search for transits of this planet 
by the discovery team did not result in a detection, although 
gra zing transits c ould n ot be ruled out. 

IRivera et al.l (120051) modeled their radial velocity data with 
self-consistent Newtonian four-body orbits assuming the three 
identified planets were in coplanar orbits. This modeling of the 
new data set yielded seemingly ti ght constraints on the inclina- 
tion of the system as foreseen by_|Laughlin & Chambers! d2001l) 



tion or the system as roreseen by Laughhn & Chambers (2UU1) 
and IRivera & Lissauerl (|2001|) . IRivera .et al l d2005l) found that 
the coplanar system inclination appeared to be ~ 50°, and the 
masses for planets b, c, and d were 2.3 Mj up , 0.2>Mj up , and 
7.5 M e respectively. Although no formal uncertainty in the incli- 
nation was given, inspection of their reported x 1 map (see their 
Fig 3) suggests a stan dard error of - 2°. 

IRivera etalJ (120051) also argued that the orbits of the two gas 
giants must have a small or even zero mutual inclination (i.e. 
they are coplanar) because the radial velocity fit quality of their 
model deteriorated when the inclinations of each of the planets 
were pushed away from 50°. However, quantification of the mu- 
tual inclination of the planets' orbits was not possible due to the 
incomplete characterization of their full three dimensional or- 
bits. Specifically missing was needed information on the planets' 
absolute orbital inclinati ons (i or 180 — i) a nd orbital longitudes 
of the ascending nodes. IRivera et al .1(120051) did not consider the 
HST astrometry or the results from lBenedict et al.l d2002l) in their 
analysis. 

Th e findings of IBenedict et al.l (|2002) and IRivera et ail 
(2005) with regards to the planets' orbital inclinations and 
masses are inconsistent, and this has led to confusion about what 
is the "best" model of the GJ 876 system. Numerous theoretical 
studies have been carried out since the discovery of the second 
planet to determine what physi cal processes gave rise to the sys- 
tem's unique architecture (e.g. ILe e & Peale 2002; Ji et al.l2002t 
iKlev et al.|[2004l 120051: IZhou et alJl2005T: ICrida et al.ll2008l) be- 
cause such specific consideration has the potential for con- 
straining general theories of planet formation and evolution. 



Determining what mechanisms could have led to the current 
system arrangement depends critically on knowing the masses 
of the planets involved, and what the current arrangement itself 
even is. Therefore, continued refinement of the GJ 876 system 
model would be valuable. 

In this paper we present new constraints on the architecture 
of the GJ 876 planetary system based on reanalysis of the previ- 
ously published high precision radial velocities and astrometry 
for the system . We aimed to reso lve t he discrepancy bet ween 
the results of IBenedict et all (120021) and lRivera et all (120051) . and 
study the degree of coplanarity in the system by the first com- 
bined analysis of the data sets using a self-consistent Newtonian 
orbit model. Our st udy represents a syne rgy in the conceptual 
advances made by Ben edict et al.l d2002l) in their analysis of 
a combined radial velocity and astrometry data set for an ex - 
oplanetary system , and t hat o f lLaughlin & Chambers! d200ll) . 
IRivera & Lissauerl d2001l) . and IRivera et al.1 (120051) iritheir use 
of Newtonian orbits to account for the signature of multi-body 
interactions in radial velocity data and infer more information 
than it is normally possible to obtain from such data. The paper 
is organized as follows. In §2 we describe the data utilized in our 
analysis. We present our analysis in §3. We conclude in §4 with 
a discussion of the implications of the results obtained and the 
potential for continued work in this area. 



2. The data 

2. 1 . Radial velocities 

We utilized the t i me ser ies radial velocities for GJ 876 presented 
by IRivera et aL I d2005l) in our analysis. These velocities were 
measured from high-resolution spectroscopic observations made 
with the HIRES instrument equipped with an iodine absorp- 
tion cell and fed by the Keck I telescope at the W. M. Keck 
Observatory. This data set contains 155 measurements that have 
a median uncertainty of 4. 1 m s _1 and were obtained over 7.6 yr. 
No adjustments were made to the uncertainties to account for 
the potential affect of stellar "jitter" (a loose term referring 
to changes in stellar spectra arising from variation of inhomo- 
geneities on the surface of star s that can be misco nstrued as a 
radial velocity change) because IRivera et al.l d2005l) achieved a 
best-fit reduced x 1 for the radial velocities very close to 1 .0, 
which indicates that there is little or no additional noise in the 
data. M ore details about these data can be found in IRivera et al.l 
(2005) and references therein. 

Other relevant radial veloci ty data for GJ 876 were pre- 
sented by iDelfosse etail (1 19981 data from t he EL ODIE and 
CORALIE spectrographs) and uvlarcv et al.1 (1200 ll data from 
Lick Observatory). We elected not to include these data sets in 
our analysis because their consideration would have been more 
confusing than illuminating. As all the radial velocities are rel- 
ative in nature, each data set included therefore requires the ad- 
dition of a free offset parameter in the orbit fitting. Furthermore, 
there is the issue of how to treat the estimated uncertainties 
when using inhomogeneous data sets. Each group has their own 
method of error estimation based on some combination of the 
photon statistics in the obtained spectra, stability of the in- 
strument used, and accuracy of the Doppler shift measurement 
method employed. Also, the uncertainties for the velocities in the 
three neglected data sets are all > 10 ms 1 and none approach 
the time baseline of the Keck velocities . Therefore, our choice 
to use only this later data set simplifies the analysis and bypasses 
a number of potential issues without ignoring very useful data. 
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We note that lRivera et al.l d2005l) also chose to focus exclusively 
on the Keck velocities for GJ 876 in their analysis. 

2.2. Astrometry 

We also included in our study the astrometry for GJ 876 that was 
obtained wit h the FGS3 instrume nt on the HST and that was 
presented by iBenedict et al.l (120021) . The data set is made up of 
observations of GJ 876 and five reference stars obtained over 27 
HST orbits and distributed in 9 epochs. The measurements span 
2.5 yr and are coincident with the Keck radial velocities. The 
majority of the observations were obtained near the periastron, 
apastron, and subsequent periastron times of planet b during one 
of its orbits (a time span of ~60d). Further observations were 
scheduled to allow breaking the degeneracy between the orbital, 
proper, and parallactic motions. 

We used the exact same data as Ben edict et al.l d2002l) with 
only one modification. We multiplied the nominal uncertain- 
ties estimated by their data reduction pipeline for the X and Y 
axis position measurements by 0.34 and 0.50 respectively. These 
re-weightings were motivated when we obtained a reduced x 1 
significantly below 1.0 for the data during preliminary model- 
ing. The weightings were iteratively adjusted to yield a reduced 
X 1 = 1.0 for the best-fit coplanar model (see below). The sep- 
arate re-weightings for the X and Y axis data are appropriate 
because the FGS instrument has two separate arms for mea- 
suring the apparent positions in the two axes. Thus, the data 
from the axes are essentially independent. The median uncer- 
tainty in the position measurements for both axes is 0.95 mas 
after re- weighting. More details about t hese data, and FGS mea- 
surements in general, can be found in IBenedict et al. I (120021) and 
references therein. 



3. Analysis 

3.1. Modeling 

Our analysis consisted of modeling the Keck radial velocities 
and HST astrometry described in §2 simultaneously. A self- 
consistent model for a four-body system (three planets and the 
host star) was used to account for the orbital motion of GJ 876 
in bot h data sets. This model was generated using the Mercury 
code (Chambers 1999) to integrate the equations of motion. All 
the bodies were assumed to be point masses and the only force 
considered was Newtonian gravity. We have previously used this 
same general method to simultaneously model radial velocities 
and eclipse times for an exoplanetary s ystem with consider ation 
of possible planet-planet interactions dBean & Seif ahrt 2008). 

We a ssumed the mass of GJ876 (m^) is 0.32 M®, as sug- 
gested bv lRivera et al. I (120051) based on consultation with empir- 
ical Mass - Luminosity relationships for low mass stars. The un- 
certainty in the estimate is probably ~10%. We did not attempt 
to account for this uncertainty because it would be prohibitively 
time consuming to repeat the analyses numerous times with dif- 
ferent assumed values. With this assumed mass, the model for 
the orbital motion of GJ 876 depended on the input masses and 
osculating orbital elements for the three planets. The orbital el- 
ements are the six usual ones: semi-major axis (a), eccentric- 
ity (e), argument of periastron (w), mean anomaly (M), incli- 
nation (/), and longitude of the ascending node (Q). The ref- 
erence epoch for the osculating e lements was taken to be HJD 
2 452 490.0 as lRivera efa l. (2005) did so that our results may be 
directly compared to theirs. Including the masses, there were 7 



parameters for each planet and 21 parameters total in the orbit 
model. 

Following [Rivera et al.l d2005l) . we always fixed the eccen- 
tricity and argument of periastron for planet d to zero. This 
planet is expected to be in a nearly circular orbit owing to tidal 
torques from the host star. Furthermore, it induces a modulation 
with semi-amplitude of only 6.5 ms -1 on the radial velocities. 
Therefore, the signature of non-zero eccentricity is negligible 
given the quality of the data and may be ignored. The treatment 
of the remaining orbit model parameters varied in the different 
analyses described below. 

Comparison of the orbit model with the radial velocities re- 
quired one additional step. A single correction factor was added 
to all the model radial velocities to shift them to the relative scale 
of the observed velocities. This offset was always a free param- 
eter in the analyses described below. The equation of condition 
for the radial velocity model was thus 

A y = RV - (y + ORBIT R ), (1) 

where RV is the measured relative radial velocities, ORBITr is 
the model radial velocity component of GJ 876's orbital motion, 
y is the offset, and A y is the residual. 

Our approach for generatin g the astrometric m odel closely 
followed the methods used bv lBenedict et al.l d2002l) . which have 
also bee n used to analyze FGS astrometry for other exoplanetary 
system s dMcArthur et al.ll2004t IBenedict et"aill2006t iBean et al.l 
12007b and binary star systems (e.g. IBenedict et all 12001). The 
model included the orbital motion of GJ 876, parallactic and 
proper motion for GJ 876 and the five reference stars, and plate 
adjustments for the 27 epochs of data. 

We selected the observations made during epoch 22 to serve 
as the astrometric constraint "plate." FGS observations of differ- 
ent stars during an epoch are carried out sequentially rather than 
simultaneously so there isn't an actual plate in the traditional 
sense. However, the sequential observations during a single HST 
orbit may be combined to form an effective plate due to stability 
of the telescope and instrument response during that time period. 
We refer to this combination of sequential observations as sim- 
ply a plate below for brevity. The choice of the constraint plate 
does not have a significant impact on th e results and our sp ecific 
choice was made for consistency with IBenedict et at] (120021) . 

The position deviations of the stars due to parallactic, proper, 
and orbital (GJ 876 only) motion were calculated in the usual 
right ascension - declination reference frame first. These devi- 
ations were then rotated about the roll angle of the FGS during 
the constraint plate observations to place them in the X - Y ref- 
erence frame of the instrument at that epoch. The equations used 



for these calculations were 

D a = P a n + fi a At + ORBIT a , (2) 

D s = P s n + fig At + ORBIT s , (3) 

= D a cos + D s sin 6», (4) 

D ri = -D a sin + D s cos 6», (5) 



where D are the motion displacements, P are the parallax fac- 
tors, n is the parallax, n are the proper motions, At is the time 
difference from the reference epoch, ORBIT are the orbital mo- 
tions, and 9 is the roll angle of the constraint plate. The a and 
6 subscripts refer to the right ascension and declination compo- 
nents respectively. The f and r\ subscripts refer to the X and Y 
components in the reference frame of the constraint plate respec- 
tively. 
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Table 1. Parameters from the coplanar analysis. 



Orbital Parameters 
Parameter 


Planet b 


Planet c 


Planet d 


m 




0.80+^^1 M ]up 


o * ^+0 i * 

8.17™;^ M e 


a (AU) 


0.20688 + "!!S 


0.13062+™™! 

v/. A^v/ui- Q 00004 


0.0208069 + °S!!n! 

u.vi.uuuuy_() .0000004 


e 


0376+ 0022 

'"-0.0019 


2657+ 00022 

u,iuJ/ -0.0017 


0.0 (Fixed) 


cj O 


184.0^ 


197.3!°;* 


0.0 (Fixed) 


M (") 


167.3^ 


311.6!J;° 


311.8+^ 


Parameter 


Value 






i (°) 


48.9+}-* 






a (°) 


251 -!e 






Fit information 
Parameter 


Value 






x 1 


860.0 






DOF 


843 






radial velocity rms (ms _1 ) 


4.2 






astrometry rms (mas) 


0.9 







Note: The orbital parameters are osculating and valid at HJD 2 452 490.0. 



The roll angle of the constraint plate was estimated to be 
26.01° + 0.05° based on comparison to ground based astrome- 
try catalogs. The actual roll angle we used for calculating the 
model was always a free parameter and the estimated value was 
treated as an observation with error to provide a constraint (i.e. 
a comparison of the estimated value to the actual value used was 
included in the overall^ 2 calculation). 

We always solved for the parallaxes and proper motions of 
GJ 876 and the five reference stars in each of the analyses de- 
scribed below. We used as observations with error the same es- 
timated parallaxes and previously measured proper m otions for 
the reference stars a lso utilized by Ben edict et alj d2002l) . Unlike 
iBenedict et al.l d2002l) . we also used the Hipparcos parallax and 
proper motion for GJ 876 (n = 212.69 + 2. 10 m as, ii a = 960.31 
+ 3.7 7 mas yr -1 , pg = -675.61 + 1.58 mas yr -1 , Perrvma n et alj 
1 19971) as observations with error. This is justified because the 
Hipparcos solution for GJ 876 is unlikely to be affected by its or- 
bital motion due to the small size (~0.3 mas) and short timescale 
(~60d) of the perturbations relative to the precision (~5mas) 
and time span (2.2 yr) of the Hipparcos observati ons of GJ 876. 

We used the same six parameter model as [Benedict et alj 
(2002) to account for changes in the plate scale, rotation, and 
offset during the different observational epochs. The equations 
of condition for the astrometry model were 

= Ax + By + C + R x (x 2 + y 2 ) ~ (f + D f ), (6) 
A, = -Bx + Ay + F + R y (x 2 + y 2 ) - (i] + £>,,), (7) 

where A are the residuals, x and y are the measured positions; 
A, B, C, F, R x , and R y are the plate parameters; and £ and rj 
nominal positions of the stars at the reference epoch. For the 
observations made during the adopted constraint epoch, the plate 
parameter A was fixed to 1, and B, C, F, R x , and R y were fixed 
to 0. These were free parameters for each of the other 26 epochs. 



The nominal positions for each star at the reference epoch were 
also free parameters. 

In total, and excluding the orbital motion component for 
GJ 876, there were always five free parameters for each of the 
six stars, six free parameters for each of 26 epochs, plus the one 
roll angle for a total of 186 astrometric only parameters. There 
were 436 FGS observations in each X and Y, as well as three 
constraints for each of the five stars and one constraint for the 
roll angle. 

We used the usual x 2 parameter as the goodness-of-fit metric 
in our analyses. The^ 2 of the model comparisons to all the data 
was calculated from Eq. (fl}, ©, and (0 along with the compar- 
ison of the certain astrometric parameters mentioned above to 
their input values. 



3.2. Coplanar study 

We earned out an analysis of the data assuming that the plan- 
ets were in coplanar orbits (i.e. we set i - it, - i c - id and 
Q. = Q.b = Sl c = Q<t)- We used a combination of grid search and 
local minimization algorithms to find the parameters that min- 
imized the x 1 between the model and the observed data. The 
parameter uncertainties were estimated by stepping out from 
the best-fit values of each parameter in turn while marginal- 
izing over the remaining parameters until the x 1 increased by 
1.0 from the minimum value. The identified orbital parameters 
and their uncertainties along with some fit quality statistics are 
given in Table Q] The identified astrometric parameters (paral- 
laxes, proper motions, and positions) are essentially the same as 
for the non-coplanar analysis (§3.3), so we give only the results 
from the later investigation because we consider them more ro- 
bust (see Tabled. 
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,1 

For the orbital orientations, we find i = 48.9 and Q = 

— 1.6 

251 0+ j^ . The inclination we determine for a coplanar model is 
very similar to that suggested by iRivera etaL I (120051) based on 
analysis of the radial velocities alone. It is also c ompletely con- 
sistent with the astrometric perturbation size that Be nedict et alj 
(2002) measured, but is not consistent with their estimated incli- 
nation for planet b (//, = 84° + 6°). The difference between our 
result and theirs arises from the radial velocity data. We have 
more and better quality data than was available then. In addi- 
tion, we have the benefit of hindsight that the dynamical model- 
ing is crucial for the radial velocity data, and that this modeling 
gives more constraints on the planets' orbital orientations than 
can usually be obtained. 

We find that the planet-planet perturbations are of such sig- 
nificance that the radial velocity data set is actually more sensi- 
tive than the astrometry to the inclination of the system in this 
unique case. This is due to the dependence of the interactions, 
which are visible in the radial velocity data, on the true masses 
of the two gas giant planets (planets b and c). The situation is 
illustrated in Fig.Q] where we show the best-fit;^ 2 for the radial 
velocity and astrometry data components along with the total 
for fixed inclinations. The perturbation due to planet b is clearly 
detected in the astrometry (with false alarm probability of 1.5 
x 10~ 5 ), but its inclusion does not alter the fit quality response 
so much that the best-fit inclination is significantly different be- 
cause of the steep response of the radial velocity fit quality. 

The astrometry fit quality component does trend lower for 
smaller inclinations, but this should not be interpreted to mean 
that the astrometry would favor smaller inclinations were it in- 
dependent of the radial velocities. This is because the astrometry 
data are not sufficient to uniquely determine all of the necessary 
planet orbital parameters. Therefore, they cannot be divorced 
from the constraints given by radial velocities and an "astro- 
metric only" inclination cannot be determined. For example, the 
astrometry data are better fitted for i = 30°, but only with the 
majority of the orbital parameters determined from the fit to the 
radial velocities. In the later case, the fit is very bad (Ay 2 > 100 
from the best-fit) and, thus, the orbital parameters determined 
for this inclination and used to generate the astrometric orbit are 
not physical. As the fit to the astrometry is still acceptable for 
i ~ 50° (we have on ly a 0.05 mas larger rms than the best-fit of 
iBenedict et alj2002l) . we conclude that the astrometry and radial 
velocities are not in disagreement and that the model giving the 
best-fit to the combined data set is the optimal one. 

Although the radial velocities are sensitive to the inclination 
of the system, there remains a degeneracy between i and 180° - i 
pairs that is unresolvable with that data alone. This is because the 
radial velocities are sensitive to the inclination through their de- 
pendence on the true masses of the planets and the masses would 
be the same for i and 180° - i. The inclusion of the astrometry 
resolves this degeneracy, and we find that the inclination of the 
system is near 50°rather than 130°. 



3.3. Constraining the degree of coplanarity 
3.3.1. Fitting the data 

As discussed in § 1 . IRivera et al.l J2005 ) found that the radial ve- 
locities are sensitive to differences in the orbital inclinations of 
planets b and c, and they reached a tentative conclusion that the 
orbits were likely coplanar. However, the mutual inclination an- 
gle (<L>) of two orbits depends not only on the inclinations, but 
also on the longitudes of the ascending nodes. In this case, the 
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Fig. 1. Relative change in the best-fit total x 2 (pluses), radial ve- 
locity component x 2 (triangles), and astrometry component x 2 
(circles) as a function of the coplanar model inclination. The in- 
set shows a region around the minimum for the total x 2 - 



mutual inclination of planets b's and c's orbits is given by the 
equation 



cos 0/, c = cos ib cos i c + sin z'j sin i c cos (Q./, - f2 e ). 



(8) 



Because of their lack of the constraint on the absolute orbital 
inclinations and orbital longi tudes of the ascen ding nodes that 
is offered by the astrometry, Rivera etlril (H005) were unable to 
quantify the mutual inclination of the two planets' orbits. 

In our case, the inclusion of the astrometry helps characterize 
the full three dimensional orbit of planet b. Therefore, we have a 
benchmark to measure misalignment relative to. This motivated 
us to use the combined data set and dynamical model to deter- 
mine the orbital orientations of planets b and c uniquely, and thus 
set limits on their orbital mutual inclination. 

To do this, we fit the data using an expanded version of the 
coplanar model. Instead of solving for a single inclination and 
longitude of the ascending node, we allowed planets b and c to 
have their own independent values (ij, Q.b, ic, and Q. c ). The data 
are not very sensitive to the orientation of planet d's orbit, so we 
set its inclination and longitude of the ascending node to that of 
planet b (i.e. id = ib and Q</ = Clb), which is the most massive of 
the three planets. 

For our initial analysis, we used the same grid search and 
local minimization technique as for the coplanar study. However, 
we found the^ 2 surface to be very chaotic. This hindered reliable 
identification of the minimum x 2 and error estimation because 
there were many valleys in the surface containing local minima 
that were statistically indistinguishable from each other. 

The difficulty with the previous method led us to ultimately 
use a Markov Chain Monte Carlo (MCMC) analysis to iden- 
tify the most likely parameter values and their confidence inter- 
vals for the non-coplanar model. The details and advantages of 
MCMC are described extensiv ely elsewhere (for discussion in a n 
astronomical context see e.g. iTegmarket al1l2004t lFordll2006h . 
Our implementation used 10 Markov chains of 10 7 points each. 
For each chain step, a jump in one of the parameters was consid- 
ered. If the jump resulted in a lower x 2 than the previous point 
the jump was accepted. Otherwise, the jump was accepted with 
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Table 2. Parameters from the non-coplanar analysis. 



Orbital Parameters 
Parameter 


Planet b 


Planet c 


Planet d 






m 


O A/1+0.1 1 A4 

2 -° 4 -om M -i"P 


n TO +0.05 A/f 

0. /8_ 003 M Jllp 


Q A 1 +0.78 ji,f 
8.41_ 075 M e 






a (AU) 


0.20700 + JJ™w 

— u.uuuuy 


0.13062 + !!™^ 

— U.UUUU3 


0.0208069 + « — 

— U.UUUUUU4 






e 


0363+ 00028 

u.ujuj 00026 


0.2683" 


0.0 (Fixed) 






(O O 


188.2+1' 

-4.0 


200.4+ \i 

— i.y 


0.0 (Fixed) 






M (°) 


163.1+*;° 


309.1+1;^ 


312.2+ 4 /' 






i (°) 


47 2+ 2 ' 2 


51 1 +36 
ji.i _ 3S 


47.2 (Tied) 






n o 


252.3*™ 


249A -H 


252.3 (Tied) 






Astrometric Parameters 
Star 


(circsec) 


i 

(arcsec) 


(mas) 


(ma.s yr _ ^ ) 


Us 

(mas yr' ) 


GJ876 


51.9001+°^ 


730.3929+™ , , 3 


215.5!";* 


955.7+j? 


-673.4!};| 


Ref-2 


-27.8463— 


767.6817+°;™ 


1 +0 ' 3 


5 5 +L6 


-16.0!|; 2 


Ref-3 


-226.2855+™™ 


759.8752" 


n 1+0.4 
•'•'--OA 


14 2 +2 ° 




Ref-4 


-297 4706 +0 0005 
z,y i .hi uu_ 0005 


639 3754+ 00005 


9 a+0.6 
^- J -0.6 


-39.8+11 


-43.0!}; 2 


Ref-5 


450 9640+ 00006 


635 2999+ 00004 


4. 7+0 6 


7 q+3.0 
'-3.3 


-1.2! 2 ; 2 


Ref-6 


351 6679+ 00005 


594 2625+ 00005 


i g+0.3 


-13 1+ 2 ' 4 




Fit information 
Parameter 


Value 










2 

X 


855.8 










DOF 


841 










radial velocity rms (ms -1 ) 


4.2 










astrometry rms (mas) 


0.9 











Note: Planet d's inclination and longitude of the ascending node were tied to those of planet b. The orbital parameters are osculating and valid 
at HJD 2 452 490.0. The coordinates (£ and rf) are relative positions in the reference frame of the constraint plate. To convert to the RA - DEC 

reference frame, the coordinates should be rotated about the inverse of the determined roll angle (26.07° + 04 o ). The y 2 and rms shown are for the 
best-fit model. 



a probability of exp(-Ay 2 /2). The characteristic jump sizes for 
each parameter were tuned to give a 20 - 40% acceptance rate. 

Each chain was initialized with a different combination of 
parameter values well dispersed from the region of parameter 
space thought to contain the lowest^ 2 . Initial tests indicated a 
typical correlation length of ~2000 points, or roughly 10 times 
the number of parameters. Therefore, we elected to record the 
parameters only every 2000 steps to save memory. Each chain 
took 30 CPU days computation time on an average desktop com- 
puter. 

All the chains converged to the same region of parameter 
space, or "burned in", within ~18 000 steps. To provide a sta- 
tistical check that the probability distributions had been th or- 
oughly sampled, we computed the iGelman & Rubinl ( 1 19921) R 
statistic for the parameter values among the chains. The statistic 
was within 10% of unity for all the parameters, which indicates 
the chains were likely long enough for robust inference. 

After trimming the burn-in points, we combined the data 
from the 10 chains to give parameter distributions with 49 910 
points. We adopted the medians of the MCMC distributions as 
the best estimates of the parameter values. The lcr uncertainties 



were taken to be the range of values that encompassed 68.3% 
of the parameter distributions on each side of the corresponding 
median. The results are given in Table [2] It should be noted that 
the errors from the MCMC analysis are correlated because they 
are calculated from the distributions arising from a simultaneous 
determination of all the parameters. This explains why the errors 
on our astrometric parameters are larger by facto r s of 2 - 3 than 
the uncorrelated errors given by Be nedict et alj d2002l) despite 
our achieving a similar astrometric fit quality (rms = 0.9 mas). 

Further support for the validity of the MCMC analysis is 
the fact that the orbital model formed by the adopted param- 
eter values (medians of the MCMC parameter distributions) is 
essentially the same as the one we initially identified as the best 
non-coplanar model using the grid search with local minimiza- 
tion technique (Ay 2 = 0.2). The main advantages of the MCMC 
method are that the parameter confidence limits account for the 
irregular ;f 2 surface, and allow calculation of composite parame- 
ter uncertainties when including correlations among the param- 
eters (see below). 

The MCMC parameter distributions for the masses, inclina- 
tions, and longitudes of the ascending nodes for planets b and c 
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Fig. 2. Probability distributions for the masses, inclinations, and longitudes of the ascending nodes for planets b and c from the 
MCMC analysis. The medians and two-sided 68.3% confidence limits are given by the dashed and dotted lines respectively. 
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Fig. 3. Probability distribution for the mutual inclination of the 
orbits of planets b and c computed from the probability distribu- 
tions for if,, Oft, z'c, and Q c using Eq.[8] The median and two-sided 
68.3% confidence limits are given by the dashed and dotted lines 
respectively. 



are shown in Fig. [2] As these plots demonstrate, the combined 
analysis of the radial velocity and astrometry data with the dy- 
namical model yielded tight constraints on the masses and or- 
bital orientations of the two planets. Using these probability dis- 
tributions, we may calculate the probability distribution for the 
mutual inclination angle between their orbits directly. The result 

is shown in Fig. [3] and we find <S>i, c = 5.0° ^ <>• 



3.3.2. Secular behavior 

The orbital parameters we identified are only valid for the refer- 
ence epoch (HJD 2 452 490.0) because the system configuration 
is varying with time. Therefore, the mutual inclination angle we 
determined is only a snapshot of the system and could be mis- 
leading about its normal characteristics. For example, we might 
have caught the system when planets b's and c's orbits were near 
their minimum or maximum mutual inclination. 

To study the time-dependency of the configuration implied 
by our model for the GJ 876 system, we integrated the planets' 
orbital motion forward for 1 Myr. A short segment of the results 
for the inclinations, longitudes of the ascending nodes, and mu- 
tual inclination for the orbits of planets b and c are shown in 
Fig. |4] We find that the orbital projection angles for the plan- 
ets are varying regularly, but with a variety of different frequen- 
cies and amplitudes. The mutual inclination between the planets' 
orbits varies with a main period of 4.8 yr and amplitude 0.15°. 
There are also two other coherent lower-amplitude periodicities 
around 60 d, which is similar to the outer planet's orbital period. 
These low-amplitude variations are slightly out of phase with 
one another and their interference leads to beat patterns over 
10 yr timescales. 

Aside from mutual inclination changes, the two planets' or- 
bital orientation angles relative to the plane of the sky librate 
with a period of 101 .9 yr. This is a projection effect arising from 
the libration of the planets' orbital nodes. Their mutual incli- 
nation is not affected by this variation and so the planets' or- 
bital nodes must be librating together. We conclude from this 
exploratory investigation that our measured mutual inclination 
for the orbits of planets b and c is likely representative of the 
long-term status of the system as the measurement uncertainties 
are an order of magnitude larger than the potential variations. 
A more thorough examination of the dynamical qualities of our 
model would be interesting, but is beyond the scope of the cur- 
rent paper. 
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Fig. 4. Long-term evolution of the projected orbital elements and mutual inclination for planets b (solid line) and c (dashed line). 
The variation is regular and continues in a similar fashion for at least 1 Myr. The bars on the left indicate the uncertainties in the 
osculating orbital elements used to initialize the simulation. 



4. Discussion 

Our analysis has revealed the full three-dimensional orbits, and 
thus the degree of coplanarity, of two planets in an exoplanetary 
system around a normal star for the first tim43- Broadly speaking, 
we find the orbits of GJ 876b and c to be coplanar. However, 
our results also imply a small, but potentially significant (~95% 
confidence), non-zero orbital mutual inclination that could be 
important. 

Dynamical friction from planetesimals during the late stages 
of planet formation is expected to result in orbits for gas gi- 
ants that are coplanar dKokubo & Idalll995l:|PoiTack et alJll99(5t 
iGoldreich et alj [2004b . Therefore, our general result provides 
further evidence for planet formation in a circumstellar disk, 
and suggests that the evolution of planetary systems might 
not lead to excitation of extreme inclinations for planetary or- 
bits relative to the original plane of the disk. Additional ev- 
idence from observations of exoplanetary systems for planet 
formation in a disk and little inclination from the original 
plane includes the coplanarity of a n exoplanet's orbit with 
a debris disk dBenedict et al.l 12006) and the stellar spin - 
planet orbit alignment o f a number of transiting planet sys- 
tems (IWinn et alJ l2005l 120061: IWolf et alJ l2007t twinn et al l 



2007; Na rita et alj|2007t iBouchv et al.ll2008t IWinn et alJl20Q: 
Loeillet et al l 120081: Uohnson et all 120081: ICochran et al.l [20081 
although see Hebrardetal. (200§for one possible exception to 
this trend. 

Our more subtle finding of a possible small degree of non- 
coplanarity in the GJ 876 system is a complement to the observa- 
tions that planetary system evolution often leads to eccentricity 
excitation and displacement of planets from their birthplace. As 
has been previously noted, the eccentricity of planet c is signifi- 
cantly non-zero (0.27), and i t is unlikely that bot h planets b and 
c could have formed in situ dLaughlin et al.ll2005l) . Therefore, it 



1 A previous c oplanarity measurement was obtained for two planets 
orbiting a pulsar (Konacki & Wolszczan 2003). 



seems probable that the system has undergone some significant 
evolutionary changes. 

iLee & Peald ([2002) have suggested that convergent migra- 
tion of planets b and c due to disk torques led to resonance cap- 
ture and eccentricity excitation. This seems to be the most likely 
explanation for the system's configuration, but there is still an 
open question of how the planets' eccentricities were kept from 
being e xcited to even higher values while the plane ts were mi- 
grating dKlev et al.ll2004ll2005tlLaughlin et al.ll2005l) . Either the 
planets actually didn't migrate very far, or there was effective 
eccentricity damping during the migration. 

Along this same line, iThommes & Lissauerl d2003l) have 
shown that resonance capture of two planets can also result in 
an inclination-type mean motion resonance that quickly leads to 
excitation of mutual inclinations of 30°. Mutual inclinations of 
60° or more can be achieved if the system experiences this si- 
multaneously with an eccentricity-type mean motion resonance. 
However, entry into the inclination-type mean motion resonance 
requires the eccentricity of the inner planet to be !~0.6, which is 
a condition that was likely not met in the GJ 876 system. Our 
finding of only a small mutual inclination for planets b and c is 
therefore a further constraint on the system's evolutionary his- 
tory. The nearly coplanar configuration of the planets' orbits 
is fully consistent with the scenario that they did not experi- 
ence the inclination-type mean motion resonance because of the 
only moderate eccentricity excitation during migration. Thus, 
the question of why the planets' eccentricities were not excited 
to higher values becomes more important. It would be interest- 
ing to investigate whether hydrodynamic simulations of differ- 
ential jmgration and resonance capture due to disk interactions 
(e.g. lKlev et al.ll2004l) could reproduce the small degree of non- 
coplanarity we have found when extended to three dimensions. 

As the GJ 876 system is the only planetary system other 
than the Solar System for which we have tight constraints 
on the degree of coplanarity it is interesting to compare the 
two. Surprisingly, we find that they share some similarities de- 
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spite their obvious differences. GJ 876b and c have a mass ra- 
tio (mb/m c = 3.38) very similar to the Jupiter - Saturn pair 
(ntjup/msat = 3.34). This seems like a coincidence because it 
is unclear how such a property of neighboring gas giants could 
be maintained in different formation environments. Gas giants 
are thought to form via runaway gas accretion on to a solid core 
so such a property would require exact timing uniquely for each 
case. 

More interestingly, Tsiganis et al. (2005) hypothesized that 
Jupiter and Saturn experienced an encounter with the 2:1 mean 
motion resonance due to migration. They suggested that this 
encounter led to eccentricity and mutual inclination excitation 
for the planets in the outer part of the Solar System, and is 
the reason for their currently non-circular and no n-coplanar or- 
bits. In contrast to the GJ 876 system though, the Tsigani set alj 
(2005) model for the Solar System has Jupiter and Saturn pass- 
ing through the resonance owing to their diverging migration. 
As a result of being caught in the resonance, GJ 876c 's eccen- 
tricity was pumped up to at least three times the value that any 
of the Solar System giant planets' orbits reach. Additionally, 
our results indicate that the GJ 876 b-c orbital mutual inclina- 
tion is potentially a few times larger as well. Furthermore, the 
Jupiter - Saturn 2:1 resonance encounter interactions also in- 
volved Uranus and Neptune, and the planets' final orbits de- 
pended on the details of the complex four body scattering. The 
GJ 876 system is known to harbor an additional low-mass planet 
in a short period orbit. It is unclear how this object was involved 
in the dynamical evolution of the system, to say nothing of other 
still-to-be-discovered planets that could potentially exist in the 
system. 

Ultimately, our determination of the orbital mutual inclina- 
tion for GJ 876b and c is just one more piece of the planet forma- 
tion and evolution puzzle. It would be useful to obtain such mea- 
surements for many other exoplanetary systems to see if most or 
all systems tend to be fairly coplanar, but with a small amount 
of mutual inclination. Systems with planets in low-order reso- 
nances are particularly interesting targets due to the constraints 
on disk interactions knowledge of their architecture provides. 
In this regard, our analysis illustrates how such measurements 
could be achieved using data obtained from existing facilities. 
The GJ 876 system is unique for the size and timescale of the 
planet-planet interactions, and so radial velocity data for other 
systems will not be as sensitive to their architecture. However, a 
number of other moderately interacting multi-planet systems are 
better astrometric targets than GJ 876 because one of the planets 
in them induces a host star perturbation larger than 1 mas, which 
is the typical measurement uncertainty of the HST FGS. If ro- 
bust astrometric characterization of one planet in a moderately 
interacting system can be obtained, then dynamical considera- 
tions could be used to constrain the degree of coplanarity for the 
other planets. 
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